Hybrid Numerical Solution of the 
Chemical Master Equation 
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ABSTRACT 

We present a numerical approximation technique for the 
analysis of continuous-time Markov chains that describe net- 
works of biochemical reactions and play an important role in 
the stochastic modeling of biological systems. Our approach 
is based on the construction of a stochastic hybrid model 
in which certain discrete random variables of the original 
Markov chain are approximated by continuous deterministic 
variables. We compute the solution of the stochastic hybrid 
model using a numerical algorithm that discretizes time and 
in each step performs a mutual update of the transient prob- 
ability distribution of the discrete stochastic variables and 
the values of the continuous deterministic variables. We im- 
plemented the algorithm and we demonstrate its usefulness 
and efhciency on several case studies from systems biology. 
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1. INTRODUCTION 

A common dynamical model in systems biology is a sys- 
tem of ordinary differential equations (ODEs) that describes 
the time evolution of the concentrations of certain proteins 
in a biological compartment. This macroscopic model is 
based on the theory of chemical kinetics and assumes that 
the concentrations of chemical species in a well-stirred sys- 
tem change deterministically and continuously in time. It 
provides an appropriate description of a chemically reacting 
system as long as the numbers of molecules of the chemical 
species are large. However, in living cells the chemical popu- 
lations can be low (e.g. a single DNA molecule, tens or a few 
hundreds of RNA or protein molecules). In this case the un- 
derlying assumptions of the ODE approach are violated and 
a more detailed model is necessary, which takes into account 
the inherently discrete and stochastic nature of chemical re- 
actions [241 EQI El [271 [M]. The theory of stochastic chemi- 
cal kinetics provides an appropriate description by means of 
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a discrete-state Markov process, that is, a continuous-time 
Markov chain (CTMC) that represents the chemical popu- 
lations as random variables [51 [TD]. If n is the number of 
different types of molecules, then we describe the state of 
the system at a certain time instant by an n-dimensional 
random vector whose i-ih entry represents the number of 
molecules of type i. In the thermodynamic limit (when the 
number of molecules and the volume of the system approach 
infinity) the Markov model and the macroscopic ODE de- 
scription are equal [^. Therefore, the ODE approach can 
be used to approximate the CTMC only if all populations 
are large. 

The evolution of the CTMC is given by a system of lin- 
ear ordinary differential equations, known as the chemical 
master equation (CME). A single equation in the CME de- 
scribes the time derivative of the probability of a certain 
state at all times f > 0. Thus, the solution of the CME 
is the probability distribution over all states of the CTMC 
at a particular time t, that is, the transient state probabil- 
ities at time t. The solution of the CME can then be used 
to derive measures of interest such as the distribution of 
switching delays ^3\, the distribution of the time of DNA 
replication initiation at different origins [^Hj, or the distri- 
bution of gene expression products [31]. Moreover, many 
parameter estimation methods require the computation of 
the posterior distribution because means and variances do 
not provide enough information to calibrate parameters [15] . 

The more detailed description of chemical reactions using 
a CTMC comes at a price of significantly increased com- 
putational complexity because the underlying state space is 
usually very large or even infinite. Therefore, Monte Carlo 
simulation is in widespread use, because it allows to gener- 
ate random trajectories of the model while requiring only 
little memory. Estimates of the measures of interest can 
be derived once the number of trajectories is large enough 
to achieve the desired statistical accuracy. However, the 
main drawback of simulative solution techniques is that a 
large number of trajectories is necessary to obtain reliable 
results. For instance, in order to halve the confidence inter- 
val of an estimate, four times more trajectories have to be 
generated. Consequently, often stochastic simulation is only 
feasible with a very low level of confidence in the accuracy 
of the results. 



Recently, efficient numerical algorithms have been developed 
to compute an approximation of the CME [16U25l [3l[5l[6l [T3l 



1321 [T] [18] . Many of them are based on the idea of restricting 
the analysis of the model during a certain time interval to 
a subset of states that have "significant" probability. While 
some of these methods rely on an a priori estimation of the 
geometric bounds of the significant subset [161 1251 [3] , oth- 
ers are based on a conversion to discrete time and they de- 
cide dynamically which states to consider at a certain time 
step [I|[i[32]. 

If the system under consideration contains large popula- 
tions, then the numerical algorithms mentioned above per- 
form poorly. The reason is that the random variables that 
represent large populations have a large variance. Thus, a 
large number of states have a significant probability, which 
renders the numerical approximation of the distribution com- 
putationally expensive or infeasible. 

In this paper we use a stochastic hybrid approach to effi- 
ciently approximate the solution of systems containing both 
small and large populations. More precisely, we maintain 
the discrete stochastic representation for small populations, 
but at the same time we exploit the small relative variance 
of large populations and represent them by continuous de- 
terministic variables. Since population sizes change over 
time we decide dynamically ("on-the-fiy") whether we rep- 
resent a population by a continuous deterministic variable 
or keep the discrete stochastic representation. Our criterion 
for changing from a discrete to a continuous treatment of a 
variable and vice versa is based on a population threshold. 

For the solution of the stochastic hybrid model, we propose a 
numerical approximation method that discretizes time and 
performs a mutual update of the distributions of the dis- 
crete stochastic variables and the values of the continuous 
deterministic variables. Hence, we compute the solution 
of a CME with a reduced dimension as well as the solu- 
tion of a system of (non-linear) ordinary differential equa- 
tions. The former describes the distribution of the discrete 
stochastic variables and the latter the values of the con- 
tinuous deterministic variables, and the two descriptions 
depend on each other. Assume, for instance, that a sys- 
tem has two chemical species. The two population sizes 
at time t are represented by the random variables X{t) 
and Y{t), where X{t) is large and Y{t) is small. Then, 
we consider for Y{t) all events Y{t) = y that have signifi- 
cant probability, i.e., Pr{Y{t) = y) is greater than a certain 
threshold. For X{t) we consider the conditional expectations 
E[X{t) \Y{t) = y] and assume that they change continu- 
ously and deterministically in time. We iterate over small 
time steps h > and, given the distribution for Y{t) and 
the values E[X{t) \ Y{t) = y], we compute the distribution 
oiY{t + h) and the values S[X(t -I- ft) \Y{t + h) = y]. Again, 
we restrict our computation to those values of y that have 
significant probability. 

To demonstrate the effectiveness of our approach, we have 
implemented the algorithm and applied it successfully to 
several examples from systems biology. Our most complex 
example has 6 different chemical species and 10 reactions. 
We compare our results with our earlier purely discrete stochas- 
tic approach and with the purely continuous deterministic 
approach in terms of running times and accuracy. 



Related Work. Different hybrid approaches have been 
proposed in the literature [121 1311 129| . As opposed to our 
approach, they focus on Monte Carlo simulation and con- 
sider the problem of multiple time scales. They do not use 
deterministic variables but try to reduce the computational 
complexity of generating a trajectory of the model by ap- 
proximating the number of reactions during a certain time 
step. The closest work to ours is the hybrid approach pro- 
posed by Hellander and Lotstedt [14]. They approximate 
large populations by normally distributed random variables 
with a small variance and use Monte Carlo simulation to 
statistically estimate the probability distribution of the re- 
maining populations with small sizes. They consider a single 
ODE to approximate the expected sizes of the large popu- 
lations. As opposed to that, here we consider a set of ODEs 
to approximate the expected sizes of the large populations 
conditioned on the small populations. This allows us to 
track the dependencies between the different populations 
more acurately. Moreover, instead of a statistical estima- 
tion of probabilities, we provide a direct numerical method 
to solve the stochastic hybrid model. The direct numerical 
method that we use for the computation of the probability 
distributions of the stochastic variables has shown to be su- 
perior to Monte Carlo simulation [5]. Another difference is 
that the method in [14] does not allow a dynamic switching 
between stochastic and deterministic treatment of variables. 

Finally, our approach is related to the stochastic hybrid 
models considered in [4] [2] and to fiuid stochastic Petri 
nets [17] . These approaches differ from our approach in that 
they use probability distributions for the different values a 
continuous variable can take. In our setting, at a fixed point 
in time we only consider the conditional expectations of the 
continuous variables, which is based on the assumption that 
the respective populations are large and their relative vari- 
ance is small. This allows us to provide an efficient numeri- 
cal approximation algorithm that can be applied to systems 
with large state spaces. The stochastic hybrid models in [4] 
[2] [IT] cannot be solved numerically except in the case of 
small state spaces. 



2. DISCRETE-STATE STOCHASTIC 
MODEL 

According to Gillespie's theory of stochastic chemical ki- 
netics, a well-stirred mixture of n molecular species in a 
volume with fixed size and fixed temperature can be repre- 
sented as a continuous-time Markov chain (X(t),f > 0) [9]. 
The random vector X(f) = {Xi{t), . . . , Xn{t)) describes the 
chemical populations at time t, i.e., Xi(t) is the number of 
molecules of type i £ {1, . . . , n}. Thus, the state space of 
X is Z" = {0, 1, . . .}". The state changes of X are trig- 
gered by the occurrences of chemical reactions, which come 
in m different types. For j £ {1, . . . , m} let Uj £ Z" be the 
change vector of the j-th reaction type, that is, u^ = u~ +u+ 
where u~ contains only non-positive entries that specify how 
many molecules of each species are consumed (reactants) if 
an instance of the reaction occurs and vector u^ contains 
only non-negative entries that specify how many molecules 
of each species are produced {products). Thus, if X(i) — x 
for some x £ Z" with x + u" being non-negative, then 
X(f + dt) = X -|- Uj is the state of the system after the oc- 
currence of the j-th reaction within the infinitesimal time 
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Figure 1: Illustration of the exclusive switch in Ex. 1 
(picture is adapted from [22]). The stochastic hybrid 
model with only three discrete stochastic states and 
two differential equations per state. 



interval [t,t + dt). 

As rigorously derived by Gillespie [10) . each reaction type 
has an associated propensity function, denoted by ai , . . . , Om , 
which is such that aj(x) ■ dt is the probability that, given 
X(f) = X, one instance of the j-th reaction occurs within 
\t, t + dt). The value aj(x) is proportional to the number of 
distinct reactant combinations in state x. More precisely, if 
X = {xi, . . . , Xn) is a state for which x + u^ is nonnegative 
then 



OiW 



Xi ■ Xt 

m = 



■(^.-1) 



if u; = (0, . 

if Uj = — Ci 
if uj = — Ci 
if u7 = -2 • 



,0), 



(1) 



ei, 



where i 7^ ^, c^ > is a constant, and e^ is the vector with 
the i-th entry 1 and all other entries 0. We set a?j(x) = 
whenever the vector x + u~ contains negative entries, that 
is, when not enough reactant molecules are available. The 
constant Cj refers to the probability that a randomly selected 
pair of reactants collides and undergoes the j'-th chemical re- 
action. Thus, if Af is the volume (in liters) times Avogadro's 
number, then Cj 

• scales inversely with A'^ in the case of two reactants, 

• is independent of N in the case of a single reactant, 

• is proportional to A' in the case of no reactants. 

Since reactions of higher order (requiring more than two 
reactants) are usually the result of several successive lower 
order reactions, we do not consider the case of more than 
two reactants. 

Example 1. We consider a gene regulatory network, called 
the exclusive switch J20j . It consists of two genes with a 
common promotor region. Each of the two gene products 
Pi and P2 inhibits the expression of the other product if a 
molecule is bound to the promotor region. More precisely, 
if the promotor region is free, molecules of both types P\ 
and P2 are produced. If a molecule of type Pi (P2) is bound 
to the promotor region, only molecules of type P\ (P2) are 
produced, respectively. We illustrate the system m Fig. Q] 
The system has five chemical species of which two have an 
infinite range, namely Pi and P2. If x — {xi, . . . ,X5) is 
the current state, then the first two entries represent the 
populations of Pi and P2, respectively. The entry x 3 denotes 
the number of unbound DNA molecules which is either zero 
or one. The entry xn (x^) is one of a molecule of type P\ 
(P2) is bound to the promotor region and zero otherwise. 



The chemical reactions are as follows. Let j G {1, 2}. 

• We describe production of Pj by DNA — >■ DNA+Pj. Thus, 



and aj (x) = Cj ■ x-j, . 



with u 



■J+2 



• We describe degradation of Pj by Pj 
and Qj+2(x) = Cj+2 ■ Xj. 

• We model the binding of Pj to the promotor by DNA + 
Pj — >■ DNA.Pj withuj+4 = — Oj— e3+ej+3 and aj+4,{x) = 

Cj+4 ■ Xj ■ X3 . 

• For unbinding of Pj we use DNA.Pj — >■ DNA + Pj with 
Uj+6{x) = Bj + 63 — ej+3 and aj+e{x) = Cj+e ■ Xj+s- 

• Finally, we have production of Pj if a molecule of type Pj 
is bound to the promotor, i.e., DNA.Pj — >■ DNA.Pj + Pj 
with Uj+s{x) = Bj and aj+s{x) = Cj+g ■ Sj+s. 



Depending on the chosen parameters, the probability distri- 
bution of the exclusive switch is bistable, i.e. most of the 
probability mass concentrates on two distinct regions in the 
state space. In particular, if binding to the promotor is likely, 
then these two regions correspond to the two configurations 
where either the production of Pi or the production of P2 is 
inhibited. 



The Chemical Master Equation. For x G Z" and i > 0, 
let p(x, t) denote the probability that the current population 
vector is x, i.e., p{x,t) = Pr{'K{t) = x). Let p(t) be the row 
vector with entries p(x,i). Given uj~, ..., u^, u^,...,u^, 
«!,..., am, and some initial distribution p(0), the Markov 
chain X is uniquely specified if the propensity functions are 
of the form in Eq. ^. The evolution of X is given by the 
chemical master equation (CME), which equates the change 
Jjp(x, t) of the probability in state x and the sum over all 
reactions of the "inflow" Qj(x — Uj)-p(x — Uj, t) and "outflow" 
aj(x) ■p(x,f) of probability [20]. Thus, 

J m 

Since the CME is linear it can be written as ^p(i) = p{t)-Q, 
where Q is the generator matrix of X with Q(x,x + Uj) = 
q:j(x) and Q(x,x) = — X^^^i c^jC^)- K Q is bounded, then 
Eq ^ has the general solution 



p(t) = p(0) ■ e'^\ 



(3) 



(Qt)' 



where the matrix exponential is defined as e*^' = ^° 
If the state space is infinite, then we can only compute ap- 
proximations of p(t) and even if Q is finite, the size of the 
matrix Q is often large because it grows exponentially with 
the number of state variables. Moreover, even if Q is sparse, 
as it usually is because the number of reaction types is small 
compared to the number of states, standard numerical so- 
lution techniques for systems of first-order linear equations 
of the form of Eq. ((21), such as uniformization [T^, approx- 
imations in the Krylov subspace [28], or numerical integra- 
tion [33], are infeasible. The reason is that the number of 
nonzero entries in Q often exceeds the available memory ca- 
pacity for systems of realistic size. If the populations of 
all species remain small (at most a few hundreds) then the 
solution of the CME can be efficiently approximated using 
projection methods |16l 1251 13] or fast uniformization meth- 
ods [5] [6] [32]. The idea of these methods is to avoid an 
exhaustive state space exploration and, depending on a cer- 
tain time interval, restrict the analysis of the system to a 



subset of states. 

Fast Solution of the Discrete Stochastic Model. Here, 

we present a method similar to our previous work [6] that 
efficiently approximates the solution of the CME if the chem- 
ical populations remain small. We use it in Section|4]to solve 
the discrete part of the stochastic hybrid model. 

The algorithm, called fast RK4, is based on the numerical 
integration of Eq. using an explicit fourth-order Runge- 
Kutta method. The main idea is to integrate only those 
differential equations in Eq. ^ that correspond to states 
with "significant probability". This reduces the computa- 
tional effort significantly since in each iteration step only a 
comparatively small subset of states is considered. We dy- 
namically decide which states to drop/add based on a fixed 
probability threshold S > 0. Due to the regular structure of 
the Markov model the approximation error of the algorithm 
remains small since probability mass is usually concentrated 
at certain parts of the state space. The farther away a state 
is from such a "significant set" the smaller is its probability. 
Thus, the total error of the approximation remains small. 
Unless otherwise specified, in our experiments we fix S to 
10~^^. This value that has been shown to lead to accurate 
approximations [6]. 

The standard explicit fourth-order Runge-Kutta method ap- 
plied to Eq. ((2| yields the iteration step [33] 

p(t + h) = p(i) + /i ■ (ki + 2 • ka -h 2 ■ ka + k4)/6, (4) 

where /i > is the time step of the method and the vectors 
ki,k2,k3,k4 are given by 



ki = p(t) • Q, 
ka = (p(t) + h ■ 



k3 = (p(i)-f/i-^)-Q, 
Q, k4 = (p(t) + h-k3)-Q. 



(5) 



Note that the entries fci(x), . . . , fc4(x) of state x in the vec- 
tors ki , . . . , k4 are given by 

m 

ki (x) = E (^i (x- Uj ) ■p(x- Uj , t) - Oj (x) ■p(x, f )) , 

fc,+i(x) = E (Qi(x-Uj)-(p(x-Uj,i)-|-/i-fc,(x-Uj)/2) 

-Q,(x)-(p(x,t) + /i-fc4x)/2)) for ie {1,2}, (6) 
fc4(x) = E (aj(x-Uj)-(p(x-Uj,i)-f/i-A:3(x-Uj)) 
-aj(x)-(p(x,f) + /i-fc3(x))). 

In order to avoid the explicit construction of Q and in order 
to work with a dynamic set Sig of significant states that 
changes in each step, we use for a state x a data structure 
with the following components: 

• a field :x..prob for the current probability of x, 

• fields x.fci, . . . , x.fc4 for the four terms in the equation of 
state X in the system of Eq. ((5|, 

• for all j with x -|- u~ > a pointer to the successor state 
X + Uj as well as the rate Oj (x) . 

We start at time f = and initialize the set Sig as the set of 
all states that have initially a probability greater than S, i.e. 
Sig := {x|p(x, 0) > 5}. We perform a step of the iteration 
in Eq. Q by traversing the set Sig five times. In the first 
four rounds we compute ki , . . . , k4 and in the final round 



1 


:hoose step size h; 


2 


'or i = 1, 2, 3, 4 do //traverse Sig four times 


3 


//decide which fields from state data structure 


4 


//are needed for ki 


5 


switch i 


6 


case j = 1: coeff := 1; field := prob; 


7 


case i G {2,3}: coeff := h/2; field := fci_i; 


8 


case i = 4: coeff := h; field := fci-i; 


9 


A..n(j . — J\..hj\, 


10 


for all X G Sig do 


11 


for j — I, . . . ,m with x + Uj > do 


12 


x.fci := x.fci — coeff ■ :>i.. field ■ aj(x); 


13 


if x+Uj ^ Sig then 


14 


Sig :— Sig U {x + Uj}; 


15 


(x-l-Uj).fci := (x-hUj).fci -1- coeff -x.field-aj (x.); 


16 


for all X G Sig do 


17 


x.prob :=x.pro6 + /i-(x.fci+2-x.fc2+2-x.fc3-|-x.fc4)/6; 


18 


x.fci :— 0; x.fca :— 0; x.fcs := 0; x.fc4 :— 0; 


19 


if x.prob < 5 then 


20 


Sig ■- Sig \ {x}; 



Table 1: A single iteration step of the fast RK4 algo- 
rithm, which approximates the solution of the CME. 



we accumulate the summands. While processing state x in 
round i, i < 5, for each reaction j, we transfer probability 
mass from state x to its successor x -f Uj , by subtracting 
a term from fci(x) (see Eq. ((6])) and adding the same term 
to fci(x -I- Uj). A single iteration step is illustrated in pseu- 
docode in Table [T] In line 20, we ensure that Sig does not 
contain states with a probability less than S. We choose 
the step size h in line 1 as suggested in (35] • In line 2-15 
we compute the values fci (x) , . . . , fc4 (x) for all x G Sig (see 
Eq. ©). The fifth round starts in line 16 and in line 17 the 
approximation of the probability p(x,t -|- h) is calculated. 
Note that the fields x.fci, . . . ,x.fc4 are initialized with zero. 

Clearly, for the solution of the CME the same ideas as above 
can be used for many other numerical integration methods. 
Here, we focus on the explicit RK4 method and do not con- 
sider more advanced numerical integration methods to keep 
our presentation simple. The focus of this paper is not on 
particular numerical methods to solve differential equations 
but rather on general strategies for the approximate solu- 
tion of the stochastic models that we consider. Moreover, 
we do not use uniformization methods as in previous work 
since uniformization is inefficient for very small time hori- 
zons. But small time steps are necessary for the solution of 
the hybrid model in order to take into account the dependen- 
cies between the stochastic and the deterministic variables. 

3. DERIVATION OF THE DETERMINISTIC 
LIMIT 

The numerical approximation presented in the previous sec- 
tion works well as long as only the main part of the prob- 
ability mass is concentrated on a small subset of the state 
space. If the system contains large populations then the 
probability mass distributes on a very large number of states 
whereas the information content is rather low since we dis- 



tinguish, for instance, the cases of having Xi{t) — 10000, 
Xi{t) = 10001, etc. In such cases no direct numerical ap- 
proximation of the CME is possible without resorting to 
Monte Carlo techniques or discarding the discreteness of the 
state space. If all populations are large the solution of X can 
be accurately approximated by considering the determinis- 
tic limit of X. Here, we shortly recall the basic steps for the 
derivation of the deterministic limit. For a detailed discus- 
sion, we refer to Kurtz [21) . 

We first define a set of functions fij such that if A'^ is large 
(recall that A*' is the volume times the Avogadro's num- 
ber) then the propensity functions can be approximated as 
q:j(x) « N ■ Pji^), where z = {z\, . . . , Zn) = x • A'"""'^ cor- 
responds to the vector of concentrations of chemical species 
and belongs to R". Recall the dependencies of Cj on the 
scaling factor N as described at the beginning of Section [S] 
For constants kj > that are independent of N , 

• Cj — kj ■ N in the case of no reactants, 

• Cj — kj in the case of a single reactant, 

• Cj — kj /N in the case of two reactants. 

From this, it follows that except for the case of bimolecu- 
lar reactions, we can construct the functions /3j such that 

a,(x) = iV./3,(z). 



A(z) = 



N 



k, 



^j' N ^ J ' ^' 



if u7 



if u, = -ei. 



(0,...,0), 



kj-Zi-zi if u- 



-Qi — ei, 



where i j^ £. In the case of bimolecular reactins (u = 
— 2 • Gi), we use the approximation 



7V./?,(z) 



kj-Xi-Zi = {\cjN)-Xi- 
Vi[xi -l) = aj{x), 



which is accurate if Xi is large, In order to derive the deter- 
ministic limit for the vector X(i) — (A'i(t), . . . , X„(t)) that 
describes the chemical populations, we first write X(i) as 

m 

X(t) = X(0) + ^u, -C.W, 

where X(0) is the initial population vector and Cj{t) de- 
notes the number of occurrences of the j-th reaction until 
time t. The process Cj{t) is a counting process with in- 
tensity aj{X.{t)) and it can be regarded as a Poisson pro- 
cess whose time-dependent intensity changes according to 
the stochastic process X(i). Now, recall that a Poisson pro- 
cess Y{t) with time-dependent intensity A(t) can be trans- 
formed into a Poisson process Y{u) with constant intensity 
one, using the simple time transform u — J^ X{s)ds, that 

is, Y{u) = y(/g A(s)ds) = Y{t). Similarly, we can describe 
Cj{t) as a Poisson process with intensity one, i.e., 



C,{t) = Y, 



a,{X{s))ds 



where Y, are independent Poisson processes with intensity 
one. Hence, for i £ {1, . . . , n} 



x^{t)^x,{o) + Y,^ 



'j» ■ ^j 



,(X(s))d. 



(7) 



where Uj = (mji, . . . , Uj„). The next step is to define Z(i) = 
X(t) ■ N~^, that is, Z(i) = {Zi{t), . . . , Zn{t)) contains the 



concentrations of the chemical species in moles per liter at 
time t. Thus, 



Zi{t) = Z,{0) + J2^i^ ■ ^~' ■ ^M / "i(^(«))'^« 



(8) 



and using the fact that aj(x) ~ N ■ /3j(z) yields 

l3,{Z{s))ds]. (9) 



Z,{t)^Z,{0) + y2uj,-N-^-Yj(N- I 
i=i ^ -^0 



By the law of large numbers, the unit Poisson process Yj 
will approach A'' ■ it at time A'^ ■ u for large A" ■ u. Thus, 
Yj {N ■ u) ^ N ■ u and hence. 



Z.{t)^Z,{0)+J2^j^- / ft(z 
.,=1 ■'o 



{s))ds. 



(10) 



The right-hand side of the above integral equation is the 
solution z(t) of the system of ODEs 



-z(t)=^u, •/3,(z(t)). 



(11) 



As shown by Kurtz [21], in the large volume limit, where the 
volume and the number of molecules approach infinity (while 
the concentrations remain constant), Z(f) — >■ z(f) in proba- 
bility for finite times i. Note that the chemical concentra- 
tions z(t) evolve continuously and deterministically in time. 
This continuous deterministic approximation is reasonable 
if all species have a small relative variance and if they are 
only weakly correlated. The reason is that only in this case 
the assumption that Zi{t) is deterministic is appropriate. 
Note that for most models this is the case if the population 
of species i is large since this implies that E[Xi{t)] is large 
whereas the occurrence of chemical reactions results only in 
a marginal relative change of the value of Xi{t). 

Example 1 (cont.). The ODEs of the exclusive switch are 
given by 

jiZl{t) = fcl ■ Z^{t) - fcg ■ 21 (i) - fcs • 21 (t) • Z3{t) 
-ffcy ■ Ziit) -f fcg ■ 24 (i) 

■^Z2{t) = k2 ■ zz{t) - ki ■ Z2{t) - ks ■ Z2{t) ■ za{t) 

+ kfi ■ 25 (i) -I- fcio ■ Z5[t) 
f^Zs{t) = -fcs • Zi{t) ■ Z-s{t) - he ■ 22(f) ■ 23(i) 

+k7 ■ 24 (i) + ks ■ 25(f) 
^24(f) = ks ■ 21(f) • 23(f) - kr ■ 24(f) 
^25(f) = ke ■ 22(f) ■ 23(f) - fcg • 25(f) 

where 21(f), 22(f), 23(f), 24(f), 25(f) denote the respective chem- 
ical concentrations. Moreover, Cj — N~^ ■ kj for j € {5, 6} 
and Cj = kj for j ^ {5, 6}. 



In [T], Ball et al. scale only a subset of the populations in 
order to approximate the behavior of the system if certain 
populations are large and others are small. Additionally, 
they take into account the different speeds of the chemical 
reactions. For a selected number of examples, they give 
analytical expressions for the distributions in the limit, i.e., 
when the scaling parameter approaches infinity. In the next 
section, we will construct a stochastic hybrid model that 
is equivalent to the one considered in [T] if we scale the 
continuous components and consider the deterministic limit. 



4. STOCHASTIC HYBRID MODEL 

A straightforward consequence of the CME is that the time 
derivative of the populations' expectations are given by 



-s[x(i)] = E'1iU.-sK(x(f))]. 



(12) 



If all reactions of the system involve at most one reactant, 
Eq. (|12|) can be simplified to 



iE[X{t)] = E"Li u, . a, {E[X{t)] 



(13) 



because the propensity functions aj are linear in x. But in 
the case of bimolecular reactions, we have either Qj(x) = 
Cj -Xi-xi for some i, I with i ^ l or Oj (x) = Cj ■ Xi ■ {xi — l)/2 
if the j-th reaction involves two reactants of type i. But this 
means that 

E[aj{X{t))]^c,-E[X,{t)-X,{t)] or 



i5[Q,(X(t))] 



E[{X,{t)f]-E[{X,im., 



respectively. In both cases new differential equations are 
necessary to describe the unknown values of E [Xi{t) ■ Xi{t)] 
and E [{Xi{t))^] . This problem repeats and leads to an infi- 
nite system of ODEs. As shown in the sequel, we can, how- 
ever, exploit Eq. (|12|l to derive a stochastic hybrid model. 

Assume we have a system where certain species have a large 
population. In that case we approximate them with contin- 
uous deterministic variables. The remaining variables are 
kept discrete stochastic. This is done because it is usually 
infeasible or at least computationally very costly to solve 
a purely stochastic model with high populations since in 
the respective dimensions the number of significant states 
is large. Therefore, we propose to switch to a hybrid model 
where the stochastic part does not contain large populations. 
In this way we can guarantee an efficient approximation of 
the solution. 

Formally, we split X(i) into small populations V(t) and large 
populations W(i), i.e. X(i) = (V(t),W(i)). Let h be the 
dimension of V(f) and n the dimension of W(i), i.e. n — 
h + h. Moreover, let D and D be the set of indices that 
correspond to the populations in V and W, respectively. 
Thus, D,D C {1, . . . , n} and \D\ = h, \D\ = h. We define 
Uj and Uj as the components of Uj that belong to D and 
D, respectively. Under the condition that V(i) — v and 
W(i) — w, we assume that for an infinitesimal time interval 
of length dt the evolution of W is given by the stochastic 
differential equation 

W(t -f dt) = W(t) -f Ej"L, Uj ■ aj (v, w) ■ dt. (14) 

The evolution of V remains unchanged, i.e., 

Pr(V(t + dt) = v + Uj |V(i) = v,W(i) = w) =aj{v,w)-dt 

The density function /i(v, w, i) of the Markov process 
{(V(i), W(i)),i > 0} can be derived in the same way as 
done by Horton et al. [17]. Here, for simplicity we consider 
only the case n = 1 which means that w = w is a scalar. 
The generalization to higher dimensions is straightforward. 
If w > then the following partial differential equation holds 
for h. 

dh{v,w,t) d{h{v,w,t) ■ Ej Uj ■Qj(v.w)) 
dt dw 

= J2aj{v-Uj,w) ■ /i(v-Uj,TO,i)-Eoj(v,™) ■ h(y,w,t). 

3 j 



\i w = Q then we have probability mass g(w, w, t) in state 
(v, lu) where 

dq(\,w,t) , , , _^ ^ 

Q^ + ^(v, w, t) ■ E, u. ■ ". (v, w) 

= Eaj(v-Uj,™) •5(v-Uj,™,t)-Eaj(v,'U') ■ giv,w,t). 
J j 

As explained in-depth by Horton et al., the above equations 
express that probability mass must be conserved, i.e. the 
change of probability mass in a "cell" with boundaries (v, w — 
dw) and (v, w + dw) equals the total mass of probability 
entering the cell minus the total mass leaving the cell. 

In order to exploit the fact that the relative variance of W is 
small, we suggest an approximative solution of the stochas- 
tic hybrid model given above. The main idea is not to com- 
pute the full density h and the mass function g but only 
the distribution of V as well as the conditional expectations 
-B[W(f) = w I V(i) = v]. Thus, in our numerical procedure 
the distribution of W is approximated by the different val- 
ues E[W{t) = w I V(t) = u], V e N" that are taken by W(i) 
with probability Pr(V{t) = v). 

Assume that at time t we have the approximation p(f) of 
the discrete stochastic model as described in Section (2] that 
is, for all states x that have a probability that is greater 
than S we have p(x, i) > and for all other states x we have 
p{'x.,t) = 0. At time t the expectations of one or more popu- 
lations reached a certain large population threshold. Thus, 
we switch to a hybrid model where the large populations 
(index set D) are represented as continuous deterministic 
variables W(f) while the small populations (index set D) 
are represented by V(i). We first compute the vector of 
conditional expectations 

*v(f) := E[W{t) = w| V(f) = v] = Ex:x=vW-p(x,f). 

Here, x is the subvector of x that corresponds to D. We 
also compute the distribution p(t) of V(f) as 

Kv,t):=Ex:x=vP(x.i)- 

Now, we integrate the system for a small time interval of 
length ft > 0. This is done in three steps as described below. 
We will write *I'v(i') for the approximation of E[W{t') = 
■w I V(t') = v]. The j-th element of the n-dimensional vector 
^v(i') is denoted by '!/'v(J, t'). The value r(v, t') denotes the 
approximation of Pr(V{t) — v) where t' £ [t,t + h). The 
vector r(f') contains the elements r(v,i'). 

(1) Update distribution. We first integrate r(t) for h 
time units according to a CME with dimension n to approx- 
imate the probabilities Pr(V{t + h) — v) by r-(v, t + h), that 
is, r(t -|- h) is the solution of the system of ODEs 

^^(^' *') - E, a, (v - u,, *._a, it')) ■ r(v - u„ i') 
-E,«.(v,*v(t'))-Kv,f') 



dt' 



with initial condition r(i). Note that this equation is as 
Eq. ([21 except that the dimensions in D are removed. More- 
over, the population sizes w are replaced by the conditional 
expectations 'l'v(t'). 

(2) Integrate. For each state v with r(v,f) > S, we com- 
pute an approximation <I>v(i+ft) of the conditional expecta- 



tion 



E[W{t + h) I V(i') = V, i' G [t, t + h)], 



that is, we assume that the system remains in state v during 
[t, t + h) and that the expected numbers of the large popula- 
tions W change deterministically and continuously in time. 
Thus, the n-dimensional vector $v(i+/i) is obtained by nu- 
merical integration of the ODE 



J.<&v(t') 



E'1iU,-a,(v,c&.(t')) 



with initial condition <l?v(i)- The above ODEs are similar to 
Eq. IHl) except that for t' G [t,t + h) the value £[Qj(X(i'))] 
is approximated by Qj(v, $v(i'))- For instance, if the j-ih 
reaction is a bimolecular reaction that involves two popula- 
tions with indices i,^ in _D then E[aj{v,W{t')) | V(f') = v] 
is approximated by Cj ■ <j)v{i, t') ■ 4'v{i, t') where the two last 
factors are the elements of the vector $v(t') corresponding 
to the i-th and ^-th population. Thus, in this case the cor- 
relations between the i-th and the ^-th populations are not 
taken into account which is reasonable if the two populations 
are large. Note that the correlations are taken into account 
when at least one population is represented as a discrete 
stochastic variable. If, for instance, i £ D and £ £ D, then 
we use the approximation Cj ■Vi-4}v{£,t') where Vi is the entry 
in vector v that represents the size of the i-th population. 

(3) Distribute. In order to approximate E\Vi/{t+h) \ V(i-f 
h)] by 'I'v(t + h) for all states v, we have to replace the 
condition V(t') = v,t' G lt,t + h) by V(i + h) = v in the 
conditional expectation $v (t+h) that was computed in step 
2. This is done by "distributing" $v(t + /i) according to the 
change in the distribution of V(t) as explained below. The 
idea is to take into account that V enters state v from v' 
during the interval [t,t + h). Assume that \t,t -\~ h) is an 
infinitesimal time interval and that q(v', v, /i), v 7^ v' is the 
probability to enter v from v' within \t,t -\- h). Then 

P(V(i+/i)=v)= E g(v',v,/i).P(V(i) = v') 

+(i-E'g(v,v',ft)).p(v(i)=v). (^^) 

Thus, we approximate E\yj[t + h) \ V(t + /i) = v] as 

E $v'(i+ft)-g(v',v,ft).P(V(i) = v'|V(t+/i)=v) (16) 

+ <i>v(i+ft)-(i-E g(v,v',;i))-p(v(t)=v|v(t+ft)=v). 

Obviously, we can make use of the current approximations 
r(i) and v{t + h) to compute the conditional probabilities 
P(V(i)=v'|V(t-|-/i) = v). For a small time step ft, <7(v',v, ft) ; 
h-aj{y' , ^v'(i)) if v' = v-f-iij and g(v', v, ft) « otherwise. 

Using Eq. (|16|) . we compute the approximation ^v(i + ft) ~ 

_B[W(i-|-ft)|V(t-f ft) = v] as 

1'v(t + ft) = 
E, <l>v-u,(t+ft) • ^^1^ ■ a,(v-u„ *._a, (i)) . ft (17) 
+*v(i+ft) ■ ^;^^(1 - E, «.(v, *v(t)) ■ ft). 
Note that the first sum runs over all direct predecessors 



Example 1 (cont.). In the exclusive switch the expected 
number of molecules of type P\ and/or P2 may becom,e high, 



K 



Figure 2: Probability distribution of the exclusive 
switch in Ex. 1. 

C(i ■'0vi(2,f) 




Figure 3: The discrete stochastic part of the stochas- 
tic hybrid model of Ex. 1. 

depending on the chosen parameters. If, for instance, a = 

C2 = Cg = Clo = 0.5, C3=Ci= C7 — Cs= 0.005, C5 = C6 = 0.01, 

and we start initially without any proteins, i.e. with prob- 
ability one in state y = (0,0,1,0,0), then after 500 time 
units most of the probability mass is located around the states 
X = (92, 2, 0, 1, 0) and x = (2, 92, 0, 0, 1) (compare the plot 
in Fig. O left). Note that xz = 0, a;4 = 1, 15 = refers to 
the case that a molecule of type Pi is bound to the pro- 
motor and xz = Xi — Q,xc, = 1 refers to the case that a 
molecule of type P2 is bound to the promotor. Since for 
these parameters the system is symmetric, the expected pop- 
ulations of Pi and P2 are identical. Assume that at a cer- 
tain time instant, both populations reach the threshold from 
which on we approximate them by continuous deterministic 
variables (we consider the unsymmetric case later, in Sec- 
tion \^. The remaining discrete model then becomes finite 
since only Pi and P2 have an infinite range in the origi- 
nal model (h = 2,h = 3). More precisely, it contains only 
3 states, namely the state vi where the promotor is free 
(xs — 1,X4, = x^ = 0^, the state V2 where Pi is bound to 
the promotor (x3, = Q,Xi = l,X5 = Q), and the state V3 where 
P2 is bound to the promotor (xs = X4 = 0,X5 = 1), see also 
Fig. QJ For i G {1, 2} let Wi{t) be the population size of Pi. 
The differential equations which are used to approximate the 
conditional expectations ^v it -f ft), j G {1, 2, 3} are 



d4>^-. {i-t') 
IP 

dp 



a — C2+i ■ (l>vi [i, t') - C4+i ■ 0V1 (*, t') 

-C2+1 ■ 0V2 (i, t') + (c? + C9) ■ (2 - i) 
-C2+i ■ (pvs {i, t') -I- (C8 4- cio) ■ (i - 1+) 



where <j)^{l,t') and (f>v{2,t') are the elements of the vector 
$vi(t'). Note that each of the 3 states vi,V2,V3 has a sys- 
tem of two differential equations, one for Pi and one for 
P2. The transition rates m the discrete stochastic part of 
the model are illustrated m Fig. [3 Thus, after solving the 
differential equations above to compute <I>vj {t -\- ft) we obtain 
the vector $v (t -|- ft) of the two conditional expectations for 
Pi and P2 from distributing $vi (i+ft), 'I'va {t+h), <I>V3 (i-f ft) 
among the 3 states as defined in Eq. (|17|l . For the parame- 
ters used m Fig. QJ left, the conditional expectations of the 
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Table 2: Results for the exclusive switch example. 



states V2 and V3 accurately predict the two stable regions 
where most of the probability mass is located. The state vi 
has small probability and its conditional expectation is lo- 
cated between the two stable regions. It is important to point 
out that, for this example, a purely deterministic solution 
cannot detect the bistability because the deterministic model 
has a single steady-state f2S)/ . Finally, we remark that in this 
example the number of states in the reduced discrete model 
is very small. If, however, populations with an infinite range 
but small expectations are present, we use the truncation de- 
scribed in Section\^to keep the number of states small. 

If at time t a population, say the z-th population, is repre- 
sented by its conditional expectations, it is possible to switch 
back to the original discrete stochastic treatment. This is 
done by adding an entry to the states v for the i-th dimen- 
sion. This entry then equals tpv(i,t). This means that at 
this point we assume that the conditional probability distri- 
bution has mass one for the value ijjv{i,t). Note that here 
switching back and forth between discrete stochastic and 
continuous deterministic representations is based on a pop- 
ulation threshold. Thus, if the expectation of a population 
oscillates we may switch back and forth in each period. 

5. EXPERIMENTAL RESULTS 

We implemented the numerical solution of the stochastic hy- 
brid model described above in C-|--|- as well as the fast solu- 
tion of the discrete stochastic model described in Section [21 
In our implementation we dynamically switch the represen- 
tation of a random variable whenever it reaches a certain 
population threshold. We ran experiments with two differ- 
ent thresholds (50 and 100) on an Intel 2.5GHz Linux work- 
station with 8GB of RAM. In this section we present 3 exam- 
ples to that we applied our algorithm, namely the exclusive 
switch, Goutsias' model, and a predator-prey model. Our 
most complex example has 6 different chemical species and 
10 reactions. We compare our results to a purely stochastic 
solution where switching is turned off as well as to a purely 
deterministic solution. For all experiments, we fixed the cut- 



3) the bistable form of the distribution is destroyed (i.e. 
promotor binding is less likely, unbinding is more likely). 

The following table lists the parameter sets (psets): 



ting threshold 5 = 10"' 
as explained in Sec. [2] 



to truncate the infinite state space 



Exclusive Switch. We chose different parameters for the 
exclusive switch in order to test whether our hybrid ap- 
proach works well if 

1) the populations of Pi and P2 are large (a) or small (b), 

2) the model is unsymmetric (e.g. Pi is produced at a higher 
rate than P2 and degrades at a slower rate than P2), 
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We chose a time horizon of t = 500 for all parameter sets. 
Note that in the case of pset 3 the probability distribution 
forms a thick line in the state space (compare the plot in 
Fig. m right). We list our results in Table [2] where the first 
column refers to the parameter set. Column 2 to 4 list the 
results of a purely stochastic solution (see Section [2]) where 
"ex. time" refers to the execution time, \Sig\ to the average 
size of the set of significant states and "error" refers to the 
amount of probability mass lost due to the truncation with 
threshold 5, i.e. 1 — X^xgSj p(^i^)- The columns 6-10 list 
the results of our stochastic hybrid approach and column 
5 lists the population threshold used for switching in the 
representations in the stochastic hybrid model. Here, "ml", 
"m2", "m3" refer to the relative error of the first three mo- 
ments of the joint probability distribution at the final time 
instant. For this, we compare the (approximate) solution of 
the hybrid model with the solution of the purely stochas- 
tic model. Since we have five species, we simply take the 
average relative error over all species. Note that even if a 
species is represented by its conditional expectations, we can 
approximate its first three moments by 

S[W(t)']^Ev(*v(i))'T(v,f) 

where the i-th power of the vectors are taken component- 
wise. Finally, in the last two columns we list the results of a 
purely deterministic solution as explained in Section [3] The 
last column refers to the average relative error of the ex- 
pected populations when we compare the purely determin- 
istic solution to the purely stochastic solution. Note that 
the deterministic solution of the exclusive switch yields an 
accurate approximation of the first moment (except for pset 
2) because of the symmetry of the model. It does, however, 
not reveal the bistability of the distribution. As opposed to 
that, the hybrid solution does show this important property. 
For pset 1 and 3, the conditional expectations of the 3 dis- 
crete states are such that two of them match exactly the two 
stable regions where most of the probability mass is located 
(see also Example 1 in Sec. [4]) . The remaining conditional 
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Table 3: Results for Goutsias' model and the predator-prey model. 
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Figure 4: Expected populations in Goutsias' model. 

expectation of the state where the promotor region is free 
has small probability and predicts a conditional expectation 
between the two stable regions. The execution time of the 
purely stochastic approach is high in the case of pset la, be- 
cause the expected populations of Pi and Pi are high. This 
yields large sizes of Sig while we iterate over time. During 
the hybrid solution, we switch when the populations reach 
the threshold and the size of Sig drops to 3. Thus, the aver- 
age number of significant states is much smaller. In the case 
of pset lb, the expected populations are small and we use 
a deterministic representation for protein populations only 
during a short time interval (at the end of the time hori- 
zon). For pset 2, the accuracy of the purely deterministic 
solution is poor because the model is no longer symmetric. 
The accuracy of the hybrid solution on the other hand is 
independent of the symmetry of the model. 

Goutsias' Model. In [11], Goutsias defines a model for 
the transcription regulation of a repressor protein in bac- 
teriophage A. This protein is responsible for maintaining 
lysogeny of the A virus in E. coli. The model involves 6 
different species and the following 10 reactions. 
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We used the following parameters that differ from the origi- 
nal parameters used in [11] in that they increases the number 
of RNA molecules (because with the original parameters, all 
populations remain small). 
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Table [3] shows the results for the Goutsias' model where 
we use the same column labels as above. We always start 
initially with 10 molecules of RNA, M, and D, as well as 2 
DNA molecules. We choose the time horizon as i = 4. Note 
that the hybrid solution as well as the purely deterministic 



solution are feasible for much longer time horizons. The 
increase of the size of the set of significant states makes the 
purely stochastic solution infeasible for longer time horizons. 
As opposed to that the memory requirements of the hybrid 
solution remain tractable. In Fig. U we plot the means of 
two of the six species obtained from the purely stochastic 
(stoch), purely deterministic (determ), and the hybrid (hyb) 
solution. Note that a purely deterministic solution yields 
very poor accuracy (relative error of the means is 95%). 

Predator Prey. We apply our algorithm to the predator 
prey model described in [5]. It involves two species A and 
B and the reactions are A — >■ 2A, A + B — ^ 2B, and B -)■ 0. 
The model shows sustainable periodic oscillations until even- 
tually the population of B reaches zero. We use this example 
to test the switching mechanism of our algorithm. We choose 
rate constants c\ = 1, C2 — 0.03, cs — 1 and start initially 
with 30 molecules of type A and 120 molecules of type B. 
For a population threshold of 50, we start with a stochas- 
tic representation of A and a deterministic representation of 
B. Then, around time 1.3 we switch to a purely stochastic 
representation since the expectation of B becomes less than 
50. Around time t = 6.1 we switch the representation of A 
because _E[j4(f)] > 50, etc. We present our detailed results 
in Table O Similar to Goutsias' model, the deterministic so- 
lution has a high relative error whereas the hybrid solution 
yields accurate results. 

6. CONCLUSION 

We presented a stochastic hybrid model for the analysis of 
networks of chemical reactions. This model is based on a 
dynamic switching between a discrete stochastic and a con- 
tinuous deterministic representation of the chemical popu- 
lations. Instead of solving the underlying partial differen- 
tial equation, we propose a fast numerical procedure that 
exploits the fact that for large populations the conditional 
expectations give appropriate approximations. Our exper- 
imental results substantiate the usefulness of the method. 
As future work we plan to include a diffusion approximation 
for populations of intermediate size. 



7. 

[1] 



[2] 



REFERENCES 

K. Ball, T. G. Kurtz, L. Popovic, and G. Rempala. 
Asymptotic analysis of multiscale approximations to 
reaction networks. The Annals of Applied Probability, 
16(4):1925-1961, 2006. 

M. L. Bujorianu and J. Lygeros. Towars a general 
theory of stochastic hybrid systems. In Stochastic 
Hybrid Systems, volume 337 of Lecture Notes in 
Control and Information Sciences, pages 3-30. 
Springer, 2006. 



[3] K. Burrage, M. Hegland, F. Macnamara, and B. Sidje. 
A Krylov-based finite state projection algorithm for 
solving the chemical master equation arising in the 
discrete modelling of biological systems. In 
Proceedings of the Markov 150th Anniversary 
Conference, pages 21-38. Boson Books, 2006. 
[4] C. Cassandras and J. Lygeros. Stochastic hybrid 
systems: Research issues and areas. In Stochastic 
Hybrid Systems, (C.G. Cassandras, and J. Lygeros, 
Ed.s), pages 1-14. Taylor and Francis, 2006. 
[5] F. Didier, T. A. Henzinger, M. Mateescu, and V. Wolf. 
Approximation of event probabilities in noisy cellular 
processes. In Proc. of CMSB, volume 5688 of LNBI, 
page 173, 2009. 
[6] F. Didier, T. A. Henzinger, M. Mateescu, and V. Wolf. 
Fast adaptive uniformization of the chemical master 
equation. In Proc. of HIBI, pages 118-127. IEEE 
Computer Society, 2009. 
[7] S. Engblom. Spectral approximation of solutions to the 
chemical master equation. Journal of Computational 
and Applied Mathematics, 229(1):208 - 221, 2009. 
[8] N. Fedoroff and W. Fontana. Small numbers of big 

molecules. Science, 297:1129-1131, 2002. 
[9] D. T. Gillespie. Exact stochastic simulation of coupled 
chemical reactions. J. Phys. Chem., 81(25):2340-2361, 
1977. 

[10] D. T. Gillespie. A rigorous derivation of the chemical 
master equation. Physica A, 188:404-425, 1992. 

[11] J. Goutsias. Quasiequilibrium approximation of fast 
reaction kinetics in stochastic biochemical systems. J. 
Chem. Phys., 122(18):184102, 2005. 

[12] E. L. Haseltine and J. B. Rawlings. Approximate 
simulation of coupled fast and slow reactions for 
chemical kinetics. J. Chem. Phys., 117:6959-6969, 
2002. 

[13] M. Hegland, C. Burden, L. Santoso, S. Macnamara, 
and H. Booth. A solver for the stochastic master 
equation applied to gene regulatory networks. J. 
Comput. Appl. Math., 205:708-724, 2007. 

[14] A. Hellander and P. Lotstedt. Hybrid method for the 
chemical master equation. Journal of Computational 
Physics, 227, 2007. 

[15] D. Henderson, R. Boys, C. Proctor, and D. Wilkinson. 
Linking systems biology models to data: a stochastic 
kinetic model of p53 oscillations. In Handbook of Appl. 
Bayesian Analysis. Oxford University Press, 2009. 

[16] T. Henzinger, M. Mateescu, and V. Wolf. Sliding 
window abstraction for infinite Markov chains. In 
Proc. CAV, volume 5643 of LNCS, pages 337-352. 
Springer, 2009. 

[17] G. Horton, V. G. Kulkarni, D. M. Nicol, and K. S. 
Trivedi. Fluid stochastic Petri nets: Theory, 
applications, and solution techniques. European 
Journal of Operational Research, 105(1):184-201, 1998. 

[18] T. Jahnke. An adaptive wavelet method for the 
chemical master equation. SIAM J. Scientific 
Computing, 31(6):4373-4394, 2010. 

[19] A. Jensen. Markoff chains as an aid in the study of 
Markoff processes. Skandinavisk Aktuarietidskrift, 
36:87-91, 1953. 

[20] N. G. V. Kampen. Stochastic Processes in Physics and 



[21 
[22; 

[23; 

[24 
[25; 

[26; 

[27 

[28; 
[29; 
[3o; 

[31 

[32; 
[33; 

[34 
[35; 



Chemistry. Elsevier, 3rd edition, 2007. 

T. Kurtz. Approximation of Population Processes. 

Society for Industrial Mathematics, 1981. 

A. Loinger, A. Lipshtat, N. Q. Balaban, and 

O. Biham. Stochastic simulations of genetic switch 
systems. Phys. Rev. E, 75(2):021904, 2007. 
H. H. McAdams and A. Arkin. Stochastic mechanisms 
in gene expression. PNAS, USA, 94:814-819, 1997. 
H. H. McAdams and A. Arkin. It's a noisy business! 
Trends in Genetics, 15(2):65-69, 1999. 

B. Munsky and M. Khammash. The finite state 
projection algorithm for the solution of the chemical 
master equation. J. Chem. Phys., 124:044144, 2006. 
P. Patel, B. Arcangioli, S. Baker, A. Bensimon, and 
N. Rhind. DNA replication origins fire stochastically 
in fission yeast. Mol. Biol. Cell, 17:308-316, 2006. 

J. Paulsson. Summing up the noise in gene networks. 
Nature, 427(6973) :415-418, 2004. 

B. Philippe and R. Sidje. Transient solutions of 
Markov processes by Krylov subspaces. In Proc. of the 
2nd International Workshop on the Numerical 
Solution of Markov Chains, pages 95-119. Kluwer 
Academic Publishers, 1995. 

J. Puchalka and A. M. Kierzek. Bridging the gap 
between stochastic and deterministic regimes in the 
kinetic simulations of the biochemical reaction 
networks. Biophysical Journal, 86(3): 1357 - 1372, 
2004. 

C. Rao, D. Wolf, and A. Arkin. Control, exploitation 
and tolerance of intracellular noise. Nature, 
420(6912):231-237, 2002. 

H. Salis and Y. Kaznessis. Accurate hybrid stochastic 
simulation of a system of coupled chemical or 
biochemical reactions. J. Chem. Phys., 122, 2005. 
R. Sidje, K. Burrage, and S. MacNamara. Inexact 
uniformization method for computing transient 
distributions of Markov chains. SIAM J. Sci. Comput., 
29(6):2562-2580, 2007. 

W. J. Stewart. Introduction to the Numerical Solution 
of Markov Chains. Princeton University Press, 1995. 
P. S. Swain, M. B. Elowitz, and E. D. Siggia. Intrinsic 
and extrinsic contributions to stochasticity in gene 
expression. PNAS, USA, 99(20):12795-12800, 2002. 
A. Warmflash and A. Dinner. Signatures of 
combinatorial regulation in intrinsic biological noise. 
PNAS, 105(45) :17262-17267, 2008. 



